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The quasiparticle spectrum of a two-dimensional d— wave superconductor in the mixed state, 
Hd <^ H <^ Hc2, is studied for large values of the "anisotropy ratio" ao = vf/va- For a square 
vortex lattice rotated by 45° from the quasiparticle anisotropy axes (and the usual choice of Franz- 
Tesanovic singular gauge transformation) we determine essential features of the band structure 
asymptotically for large an, using an effective one-dimensional model, and compare them to numer- 
ical calculations. We find that several features of the band structure decay to zero exponentially fast 
for large ao- Using a different choice of singular gauge transformation, we obtain a different band 
structure, but still find qualitative agreement between the ID and full 2D calculations. Finally, we 
distort the square lattice into a non-Bravais lattice. Both the one- and two-dimensional numerical 
calculations of the energy spectra show a gap around zero-energy, with our gauge choice, and the 
two excitation spectra agree reasonably well. 



I. INTRODUCTION 

The quasiparticle spectrum of wave superconduc- 
tors in the mixed state, Hd H ^ Hc2, has been 
the object of recent detailed studies. High-Tc supercon- 
ductors are extreme Type-II superconductors and in a 
magnetic field H ^ Hd develop a vortex lattice, the ge- 
ometry of which depends somewhat on the strength of 
the magnetic field. The d— wave symmetry implies the 
existence of four points on the zero-field Fermi surface 
where the gap vanishes. This fact has important implica- 
tions for the low-temperature thermodynamics of these 
systems, as there are states available even at very low 
temperatures. 

The main questions that have been addressed in the 
literature are the fate of these gapless points in a finite 
magnetic field and, more generally, the nature of low- 
energy states of the quasiparticle spectrum. The stan- 
dard approach to studying the quasiparticle spectrum in 
a superoanductor is through the Bogoliubov-de Gennes 
equatioru, which allows for a spatially varying order pa- 
rameter in a natural way. Gor'kav and Schrieffero and, in 
a more recent paper, AndersonB neglected in a first ap- 
proximation the spatially dependent superfluid velocity 
in the vortex lattice. They predicted that the quasipar- 
ticle spectrum in a magnetic field H <^ Hc2 is character- 
ized by broadened Landau levels. However, their key as- 
sumption, namely treating the superfluid velocity in peij=. 
turbation theory, is questionable as work by Mel'nikovEl 
shows that this singular perturbation strongly mixes Lan- 
dau levels. 

A new approach to tackle these questions was pio- 
neered by Franz and TesanovicH. They introduced a sin- 
gular gauge transformation that takes into account the 
supercurrent distribution and the magnetic field on an 
equal footing. Starting from the linearized Bogoliubov- 



de Gennes cquationB in a magnetic field, they mapped 
the original problem onto one of diagonalizing a Dirac 
Hamiltonian in an effective periodic vector and scalar po- 
tential with vanishing magnetic flux in the unit cell. The 
new Hamiltonian is more suited to numerical diagonal- 
ization, using standard band structure calculation tech- 
niques. Within their model, Franz and Tesanovic found 
that the low-energy quasiparticle spectrum stays gapless 
at the center of the Brillouin zone (F point), even in the 
presence of a magnetic field. They also noticed that, for 
large anisotropy ratio an = vf/va (where va = Aq/pf 
is the quasiparticle velocity parallel to the Fermi surface 
and Aq is the maximum magnitude of the bulk zero-fleld 
gap function), more low-energy states arise. For larger 
values oi an, they found entire lines in the Brillouin zone 
where the energies appeared to vanish, within the uncer- 
tainty of their numerical calculations. 

Because it is difficult to decide by direct numerical cal- 
culations whether there are actual zeroes of the. disper- 
sion relations, we undertook in a previous papeiu to em- 
ploy symmetry considerations, as well as a perturbative 
analysis and alternative calculational methods to eluci- 
date the properties of solutions of the Franz-Tesanovic 
(FT) equations. The present work continues these ef- 
forts by analysing a one-dimensional variant of the FT 
model, wkich was also considered by Knapp, Kallin and 
BerlinskyQ. 

The one-dimensional model has the advantage that it 
can be solved numerically or analytically with arbitrary 
accuracy, and its asymptotic behavior, for large values 
of ao can be readily extracted. Moreover, the behavior 
of the ID model and the two-dimensional FT model are 
very close to each other for large values of the anisotropy 
ratio ao- Thus we can gain valuable insight into the 
solutions of these equations. 

We note that relatively large values of the anisotropy 
ratio are relevant to actual high-Tc superconductors. 
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namely an ranging between 10 and 20. (From angle- 
resolved photoemission spectroscopy and thermal con- 
ductivity measurements, the value of a/) for high-X, su- 
perconductors, faarns out to be about 14 for YBCOo and 
20 for Bi2212yM) 

The singular gauge transformation introduced by 
Franz and Tesanovic is exact in the region of space out- 
side the vortex core where the magnitude of the gap pa- 
rameter |A(r)| is independent of position. As we are in- 
terested in ther situation of weak magnetic fields, where 
the vortex cores occupy a very small fraction of the total 
area, it is natural to suppose that errors in the energies 
introduced by an arbitrary treatment of the core interior 
would not be significant. However, the situation is mose 
subtle, as was noted by Vafek, Melikyan and Tesanovioll. 
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FIG. 1. (a) A and B sublattices and vortex lattice unit cell, 
(b) The corresponding magnetic Brillouin zone. 

The energy states of interest have a large weight in 
the vicinity of the vortex cores, and the energies can 
be shifted by amounts which are finite on the scale of 
El = hvp/d {d is the distance between nearest vortices 
belonging to the same FT sublattice, as defined in Fig. |^), 
depending on what one assumes for the behavior inside 
the vortex core. Alternatively, one can introduce a math- 
ematical boundary condition at the radius of the vortex 
core. The wavefunction, which obeys the FT equations 
outside the core, must satisfy appropriate boundary con- 
ditions at the boundary of the core. It appears, however, 
that there are a variety of possible boundary conditions 
which remain mathematically sensible and distinct in the 
limit of vanishing core radius. An analysis of the behavior 
inside the core is therefore ultimately necessary to find 
the behavior appropriate to a given microscopic Hamil- 
tonian or for an actual high- Tc superconductor, even in 
the limit of weak magnetic fields. 

We do not address this problem here. Instead, we fol- 
low the approach, of Franz and TesanovicO and our own 
subsequent workQ, in which, for a given choice of the sin- 
gular gauge transformation, the vortex core is treated by 
smoothing the superfluid velocities over a small area. Al- 
though this approach appears to be well defined in the 
limit where the core-radius vanishes relative to the vortex 



separation, the energy spectrum obtained does depend 
on the particular choice made in the FT gauge trans- 
formation, (i.e., the energy spectrum depends on which 
vortices are assigned to each of the two FT sublattices.) 
Similarly, each of the one-dimensional models considered 
in the present paper depend on the choice of FT gauge 
used in the 2D model from which it was derived, and 
the resulting dispersion relations reflect the differences 
between the parent 2D models. 

The essence of the one-dimensional model is to smear 
the vortices in the y direction (as defined in Fig. |l]). This 
direction is the "hard" direction, meaning that wavefunc- 
tions in this direction change very slowly. In particular, if 
we focus on the low-energy states, any fluctuation of the 
quasiparticle wavefunctions in the y direction (for large 
anisotropy) will be energetically very costly. For this 
reason, it is reasonable to assume that asymptotically 
for a£) — > cxD the lowest-energy states will be transla- 
tionally invariant in the y direction (and, more generally, 
they will be represented by plane waves in the y direction 
with a wavevector ky equal to a reciprocal vector of the 
vortex lattice, if higher energy states are to be taken into 
account). Of course, this assumption depends crucially 
on the orientation of the vortex lattice with respect to 
the quasiparticle anisotropy axes. We consider a square 
vortex lattice tilted by 45° with respect to the anisotropy 
axes. If the angle was different, a line going through a 
vortex along the "hard" direction would not necessarily 
intersect other vortices (incommensurate case) or could 
intersect another vortex some number of unit cells away 
from the vortex we considered (commensurate case). In 
either case the description would be more complicated 
than the one, considered here. 

Mel'niko-vQ was the first to realize that the Bogoliubov- 
de Gennes equation for a d— wave superconductor in the 
mixed state is simplified in the large anisotropy limit 
ao ^ 1. He derived a one-dimensional approximate 
model and described how to solve it, although he con- 
fined his analysis to the semiclassical vj3rsion of these 
solutions. Knapp, Kallin and Berlinksytl, starting from 
the FT equations, derived a ID model valid for large an 
and carried out a fully quantum mechanical calculation 
of the properties of this model. They worked ki the sin- 
gular gauge introduced by Franz and TesanovicH matched 
the solutions in terms of parabolic cylinder functions at 
the boundaries of the unit cell to obtain an exact ex- 
citation spectrum for the one-dimensional Hamiltonian. 
This matching, although exact, was nevertheless carried 
out numerically and therefore did not give analytical in- 
formation on the dependence of the features of the band 
structure on the anisotropy aiy. 

We have found that the most important features of 
the one-dimensional model can be understood and cal- 
culated using an approach which emphasizes the under- 
lying physics in a more transparent way. In essence, the 
low-energy states are localized around one of the two vor- 
tices in the unit cell and, in the large anisotropy limit, 
the overlap between linearly independent wavefunctions 
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can become exponentially small as a function oi an if 
the potential wells these states are confined in are spa- 
tially separated. In the singular FT gauges that exhibit 
this exponential behavior, features of the energy bands, 
such as the minimum energy in the FY direction (away 
from the F point) or the bandwidth of the FX lowest en- 
ergy band, can be calculated analytically in the WKB 
approximation. 
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(a) (b) 
FIG. 2. (a) ABAB vortex lattice unit cell (plotted with 

four vortices per unit cell), (b) A ABB vortex lattice unit 

cell. 

The lower dimensionality of the model under study in 
this paper, allows also for much smaller scale numeri- 
cal calculations to find the quasiparticle spectrum than 
the two-dimensional case, as was already pointed out 
by Knapp, Kallin and BerlinskyB. Some of the features 
which were hard to resolve in the full two-dimensional 
problem, such as the band structure along the FX sym- 
metry line and the apparent lines of zero-energy states 
that formed parallel to the FX line and intersected the 
FY line at non-symmetry points, become readily accessi- 
ble in the one-dimensional model. In this work we have 
performed such numerical calculations for two different 
singular gauge choices and found good agreement with 
the asymptotic analysis and the data available for the 
two-dimensional model in both cases. In the first sin- 
gU|La.r gauge, which we called ABAB following Vafek et 
alX^ (see Fig. ||) , the low-energy features of the spectrum 
exhibit the above-mentioned exponential dependence on 
au: while in the second gauge, which we called AABB, 
they do not. 

If, after including the appropriate vortex core physics, 
the correct energy spectrum in the continuum limit turns 
out to be related to the results in the ABAB gauge, we 
could use the above results to improve our understand- 
ing of the relation between the semiclassical description 
of the quasiparticle spectrum in the vortex state, as de- 
rived by Volovikllj and the fully quantum-mechanical pic- 
ture. The semiclassical density of states is characterized 
by two crossover energy scales, Ei and E2 {E2 < Ei). 



For energies E Ei = hvp/d the density of states is, 
on average, quantitatively identical to that in the ab- 
sence of a magnetic field, while for E < Ei the density 
of states can be evaluated in a semiclassical approxima- 
tion and turns out to be constant. The semiclassical 
description breaks down at the energy scale E2, where 
a full quantum mechanical|-calculation becomes neces- 
sary. In our previous paperQ, we noticed that the value 
of the crossoveiL energy scale E2 predicted by Kopnin 
and VolovikllXO E^^ « hvA/d did not seem to match 
with our numerical findings, at least for the vortex lat- 
tice geometries we considered. In fact, it seemed that E2 
went to zero much faster |than linearly in l/ao- Further- 
more, recentpSpecific heattZI and low temperature thermal 
conductivityliS experiments have been performed in the 
regime E < i?!"^ and still show good agreement with 
the semiclassical predictions, which seems to suggest that 
the correct crossover scale between the quantum mechan- 
ical and semiclassical regimes is much smaller than E^^ , 
for large anisotropy ratios. The one-dimensional model 
provides the key to this puzzle, as both the bandwidth 
of the lowest FX energy band and the minima of the 
FY band (away from the F point) are shown to be sup- 
pressed exponentially for large a d , instead of linearly. In 
particular, the larger of these two energies will determine 
the crossover between the semiclassical and quantum me- 
chanical regimes, where the details of the band structure 
become important to determine the density of states. 

We also studied the one-dimensional energy spectrum 
of a non-Bravais vortex lattice, with the usual choice 
of singular gauge [ABAB). In this case, we previously 
foundEI that the two-dimensional band structure becomjea 
gapped. It has been suggested recently by Vishwanathlla 
that this gap may be a spurious feature of the Franz- 
Tesanovic ABAB model; indeed one can make other 
choices for the singular gauge transformation, assigning 
the vortices to sublattices in a different way, so that the 
spectrum is gapless for the very geometry in question. 
We do not try to answer here the question of whether a 
gap would be present in this geometry when one treats 
properly the region inside the vortex cores. We limit 
ourselves here to the observation that given a particular 
choice of the singular gauge transformation the spectra 
of the two-dimensional model and of the one-dimensional 
model derived from it agree reasonably well at long values 
of In both cases the lowest positive and negative en- 
ergy bands are separated by a finite energy gap, and the 
minimum energy separation between the two bands does 
not occur at the F point (for anisotropy au 7^ !)■ Simi- 
larly to other results discussed in this paper, this energy 
gap could be a feature which depends on the choice of 
vortex core physics. Whether it survives or not the cor- 
rect choice of boundary conditions is not within the scope 
of this paper. We limit ourselves to note that given a sin- 
gular gauge transformation, we find good agreement be- 
tween the spectra of the one- and two-dimensional Hamil- 
tonians. 

To summarize the structure of the paper, the one- 
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dimensional model Hamiltonian is derived in Sec. ||. 
The asymptotic analysis of the band structure for large 
anisotropy follows in Sec. [II. Numerical results can be 
found in Sec. IV, where a comparison between two dif- 
ferent singular gauge choices is shown explicitly, together 
with a study of the excitation spectrum of a non-Bravais 
vortex lattice. Conclusions follow in Sec. 0. A deriva- 
tion of the WKB approximation and an account of the 
numerical methods can be found in Appedix A and B, 
respectively. 



II. DERIVATION OF THE ONE-DIMENSIONAL 
MODEL 

The properties of the quasiparticle spectrum in a 
superconductor in an external magnetic field are best 
studied in the framework of the Bogoliubov-de Gennes 
equation^. To establish notations and for overall clarity, 
we will very briefly summarize the most important points 
of the derivation of the Bogoliubov-de Gennes equation 
for a d-wave superconductor in the miged state, r^ore 
details con be found in our earlier paperQ (see alsal3). 

We will consider the Bogoliubov-de Gennes operator 
for a (ia;y-superconductor instead of the more conven- 
tional dx2_y2, purely for notational simplicity; results do 
not depend on this choice. Also, because we are only 
interested in the low-energy properties of the spectrum 
and because we consider magnetic fields H <C Hc2 (i.e. 
such that the size of the vortex lattice unit cell is much 
larger than 1/kp), we can linearize the Bogoliubov-de 
Gennes equation around one of the four gapless points on 
the Fermi surface. Choosing p = (0,^^), the linearized 
Bogoliubov-de Gennes operator readaS 



" I ^{pl, A*(r)} ~MPy + lAy) ) ' 

where A(r) = Aoe*'^*^''^ is the Ginzburg-Landau order 
parameter and the brackets represent symmetrization 
{a, b} = \{ah -\- ha). 

In the |-| singular gauge introduced by Franz and 
TesanovicH, we can eliminate the position dependent 
phase factor e"^*^'^) from the off-diagonal components of 
the Bogoliubov-de Gennes equation (|l|). The resulting 
Bogoliubov-de Gennes operator after the gauge trans- 
formation is 
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where wa = ^o/pf- The superfluid velocities corre- 
sponding to the A and B sublattices of the vortex lattice, 
as defined in Fig. |[ are 
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where (j}{r) — (j)A[r) + cf)B{r). The operator (g) describes 
the dynamics of a non-interacting Dirac particle in a pe- 
riodic scalar and vector potential and general featurgs-of 
its spectrum have been studied in detail elsewhereO'DO. 

The Bogoliubov-de Gennes equation should in prin- 
ciple be solved self-consistently, thus finding the order 
parameter distribution (from which the geometry of the 
minimum free-energy vortex lattice configuration can be 
inferred) together with the energies and wavefunctions of 
the quasiparticle excitations. At low temperatures and 
for the magnetic field range considered here, to a good 
approximation, the order parameter distribution can be 
computed in Ginzburg-Landau theory and substituted in 
the Bogoliubov-de Gennes equation, without iterating it 
to achieve self-consistency. We studied the quasiparticle 
spectrum in this approximation, choosing a square vortex 
lattice rotated by 45° from the quasiparticle anisotropy 
axes. 

In this work we want to focus on the large anisotropy 
ctD = vf/va S> 1 limit. For this reason we have de- 
cided to study a one-dimensional model which captures 
the essential physics of the large- anisotropy regime and 
we believe becomes asymptotically exact in the — > oo 
and low-energy limit. The essence of the model is to 
smear out the vortices in the y direction (the "hard" 
direction). The absolute value of the low-energy two- 
dimensional wavefunctions has only small ripples and 
overall fluctuations in the "hard" direction, thus motivat- 
ing the claim that asymptotically the low-energy eigen- 
states of (||) should factor into a plane wave in the y 
direction times a Bloch state in the x direction. 

To simplify equations we will write lengths in units 
of the distance between nearest-neighbor vortices d and 
energies in units of Ei = Tivp/d. Let us define the di- 
mensionless periodic potentials 



U. 



A,B 



md 



^ sx,y 



dx.'. 



ed 

he " 



(4) 



where the vector potential is chosen to satisfy the Lan- 
dau gauge Ax = —dBy and Ay = 0. The Bogoliubov-de 
Gennes operator can then be rewritten in dimcnsionless 
form 



— dx + ^ (U^ 



(5) 
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The one-dimensional approximation amounts to as- 
suming that gauge invariant quantities can only be a 
function of x. In particular, the potentials t/, 
U^'^ are gauge invariant quantities and thus 



jf^-s and 



dy (dxcj) 



A.B 



2mj) = 



which, away from the vortices, implies 
dy(f>^'^ — —2ttx + const. 



(6) 



(7) 



The constant of integration in the above expression 
changes discontinuously when crossing a line of vortices 
and can be determined by imposing that the potentials 
U:^'^ are periodic (with period 1) and have zero average 
and that 



i 



Wct)"^'^ ■ dl = 2t: Na,b, 



(8) 



where Na,b is the number of A or B vortices in the unit 
cell. The path C is a rectangle enclosing a line of vortices 
with infinitesimal width in the x direction and spanning 
the vortex lattice unit cell in the y direction. Further- 
more, in the large-anisotropy limit, vf^ and will both 
vanish. The periodicity of the potentials can be used to 
further simplify the calculation by writing the wavefunc- 
tions in Bloch form 
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The final result is that the one-dimensional effective 
Hamiltonian acting on Bloch states takes the form 
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We have considered three different vortex arrange- 
ments, two of which represent the same physical situ- 
ation of a square Bravais lattice of vortices and one is an 
inversion-symmetric non-Bravais lattice. The two Bra- 
vais vortex lattices differ in the labeling of th£_.vortices, 
as can be seen in Fig. |[ Following Vafek et alS3, we will 
refer to them as the ABAB and AABB configuration, 
respectively. 




FIG. 3. The one-dimensional periodic potentials U^{x) 
and U'^ix) in the ABAB gauge are plotted in the unit cell. 
Energies are in units of hvp/d. 



In the ABAB case, we can work with a unit cell with 
one A and one B vortex. In this case the potentials 
{x) and (x) can be calculated as outlined above to 
find 
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-2'KX — ^TT, — ^ < a; < ^ 
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(11) 



(12) 



The [/-^(x) and U^{x) potentials in the ABAB gauge 
are plotted in Fig. ^. 

Alternatively, we considered the AABB configuration, 
in which case one has to work with a unit cell contain- 
ing four vortices (two of type A and two of type B). In 
this case the discontinuity of the potentials across a line 
of vortices is tt (this should be compared to 2t: in the 
case of the ABAB lattice) and the C/^(x) and U^{x) 
potentials are identical across the unit cell 



W'ix) = U"{x) = 
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The one-dimensional model can be derived also for a 
class of non-Bravais lattices. In particular, we studied 
an inversion-symmetric lattice, whose quasiparticie spec- 
trum we analyzed in detail in a previous papeiu. This 
lattice has two vortices per unit cell whose distance along 
the diagonal, 2Rq (see Fig. ^, is different from (V2/2)d 
(which corresponds to evenly spaced A and B vortices, 
necessary condition for a Bravais lattice). For concrete- 
ness, let us consider the case where the two sublattices 
A and B are still square lattices with spacing d, but 
the distance between nearest A and B vortices is now 
{2xyy/2)d. Switching back to dimensionless units, the 
A and B vortices in the unit cell are located at coor- 
dinates {—Xtj,—Xy) and {xv,Xv) respectively, instead of 
(±1/4, ±1/4), as in the Bravais lattice case. Here we 
are particularly interested in determining whether the 
appearance of a gap in the excitation spectrum (as was 
found before in the two-dimensional numerical analysis) 
is peculiar to the two-dimensional model or is a feature 
of the one-dimensional approximation as well. To this 
end, we have calculated the U'^{x) and U^{x) potentials 
when the lines of vortices are located at a; = ±x„ 



U^ix) 
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— 27ra; — {2xy 



-2'KX + (2Xy 
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1)k, 
1)k, 



1)k, 
1)k, 



-\<x< 



^ < X < Xy 



X 1) X 



(14) 



(15) 



Note that the discontinuity of the U'^{x) and U^{x) po- 
tentials at the vortex lines is always 27r, regardless of the 
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location of the vortices. These potentials will be used 
in the analysis of the excitation spectrum, both in nu- 
merical calculations and in the asymptotic analysis of 
prominent features of the band structure for the ABAB 
vortex lattice. 



III. ASYMPTOTIC ANALYSIS OF THE BAND 
STRUCTURE FOR THE ABAB VORTEX 
LATTICE 

The one-dimensional model can be solved exactly 
by matching the parabolic cylinder function solutionscl 
across the boundary of the unit cell. This method re- 
quires a numerical exact calculation to match the bound- 
ary conditions and leads to results which are practically 
indistinguishable from the numerical diagonalization of 
the onc-dimensional model, as shown by Knapp, Kallin 
and BerlinskyO. Here we use an alternative approach 
which does not provide an analytical solution of the one- 
dimensional model for arbitrary values of the anisotropy 
but it allows us to understand in a relatively simple way 
the leading behavior of the quasiparticle spectrum for 
large values of the anisotropy ratio a^. In this section 
we will limit our analysis to the ABAB vortex lattice, al- 
though we will comment on the other gauge choice that 
was introduced in the previous section. 

In essence, in the large anisotropy limit and for the 
ABAB configuration, the low-energy wavefunctions are 
localized in the triangular potential wells U^{x) -\- ky and 
{x) — ky , defined in equation ( p^ with the potentials 
( |l^ ) . The tails of these wavefunctions overlap in the clas- 
sically forbidden regions of the potential wells, and the 
overlap becomes exponentially small for large ajj . Using 
1/aD as our "small" parameter, we can develop a WKB 
approach which can be used to calculate the leading ex- 
ponential dependence on an of the tails of the wavefunc- 
tions. This, in turn, allows us to understand the features 
of the band structure, like the lowest energy band widths 
and energy splittings at near crossings. 

There are a number of consequences of the small over- 
lap between wavefunctions localized around different vor- 
tices, as the anisotropy grows larger. First of all, note 
that the one-dimensional Hamiltonian (10) can be split 
in the following way 
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ID 



kyCT'i 
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where Tip, the one-dimensional Hamiltonian computed at 
fcj^ = 0, is independent of ky and 0-3 is the diagonal Pauli 
spin matrix. This implies that the expectation value of 
TiiD in the state ip(x), £ =< Hid >^/>, obeys 



d£ 

dk,. 



(17) 



non-zero component will be the upper (resp. lower) one. 
In this case, < (73 >^= -1-1 (resp. -1), and the slope of 
the energy bands is fixed, regardless of the value of ajj . 
This is the case for lowest energy eigenstates close to the 
Y point, especially at large anisotropy, as the potentials 
that go into the one-dimensional Hamiltonian (^o|) cor- 
respond to a roughly triangular well for one component, 
with the bottom of the well close to the energy eigen- 
value, and an almost constant potential for the other 
component, with a fairly large energy difference between 
the constant potential and the low-energy eigenvalue. In 
this approximation, the low-energy eigenstates of TiiD 
are Airy functions, for one component, and zero for the 
other. In the repeated band scheme, we expect the band 
structure close to the Y point to look like parallel lines, 
the positive energy ones with slope -1-1, and the nega- 
tive energy ones with slope -1 (the negative slope curves 
with positive energy in Fig. correspond to values of 
— 27r < ky < — TT, in the extended zone scheme). This 
is because around ky — tt, the potential U^{x) + ky is 
pushed up in energy, while U^{x) — ky is pushed down, 
and therefore, in the limit ao 00, the positive energy 
states are strongly localized around the left vortex in the 
unit cell, while the negative energy states are localized 
around the right hand one. 




If the wavefunction ip{x) is localized entirely around the 
left (respectively right) vortex in the unit cell, its only 



FIG. 4. Band structure of the one-dimensional 
in the ABAB gauge for a square lattice with anisotropy 
an = uf/wa = 12 along the FY direction. The two branches 
of each band correspond to < fcy < vr (solid curves) and 
— 2tt < ky < — TT (dashed curves). Only positive energy 
bands are plotted for clarity, negative energy bands can be ob- 
tained through particle-hole symmetry. Energies are in units 
of El — hvp/d. 

Let us turn our attention to the first near crossing at 
approximately zero energy along the FY line, closest to 
the Y point. We claim that the energy splitting between 
the first positive and negative energy bands is exponen- 
tially small in the anisotropy ao- We want to calculate 
the leading behavior of this splitting for large an , there- 
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fore we will work in the WKB approximation. The lead- 
ing exponential expression of the WKB wavefunctions 
corresponding to energy E takes the form 

^^^-^ _ ^±an P ^-{E-ky-U^{x)){E+ky-UB{x)) ^-^g^ 

and analogously for the lower component as dis- 



cussed in more detail in Appendix A. We can calculate 
the lowest energy levels at the Y point using the Bohr- 
Sommerfeld quantization rule. If \E\ < 27r, the turning 
points for the left vortex (centered at x — —1/4) are 
x^ — — J — ^ and x^ — —1/4, therefore the quantiza- 
tion condition reads 



with n > 1. On the right hand side we have used the 
connection rules for a soft boundary on one side and a 
hard boundary on the other side of the classically allowed 
region, as we assumed that the energy of the eigenstate 
is small compared to the height of the triangular well. 
The integral on the left hand side of equation (|l^) can 
be evaluated explicitly to find 



i(^E + ^)^E{ET^-^sinh~ 



1 



■.E^'\ 



(20) 

where we have used sinh~^ x ^ x for a; ^ 1. This con- 
dition translates in our case to E <^ which, for large 
enough au, is satisfied even for large n, thereby justifying 
the usage of the Bohr-Sommerfeld quantization formula. 
In conclusion, we find the quasiparticle spectrum at the 
Y point to be 

with n = 1,2,... From our previous discussion, we also 
know that as we go away from the Y point the slope of 
the energy band is either +1 or -1, depending on their 
sign. This means that, if we neglect the level repulsion 
for the moment, the lowest positive and negative energy 
bands will cross at zero energy (because of the underlying 

2/3N 



particle-hole symmetry) at k* 



TT 1 



We can now find the magnitude of the splitting be- 
tween the lowest positive and negative energy bands at 
k*. Our method follows closely the LCAO approximation 
in molecular physics. We consider two wavefunctions cor- 
responding to the lowest-energy states localized around 
each vortex for ky = k*, that is spinors with either a non- 
vanishing upper or lower component. If these states were 
true eigenstates of the one-dimensional Hamiltonian (p^), 
the corresponding eigenvalue would be zero, by definition 
of k* . Because the "small" component of the spinor wave- 
function is not exactly zero, we can evaluate variationally 
the exponential dependence of the lowest eigenvalues on 
the spatial separation of the wavefunctions calculating 
the overlap integral between the WKB wavefunctions of 
the two zero-energy states. Because of the periodicity 
of the potential, there are in principle two regions which 
are classically forbidden and therefore contribute to the 
tunneling amplitude between the wavefunction localized 
around one vortex and the other. On the other hand, 
for low-energy states, one of these regions (going from 
one vortex site to the other, regardless of the value of 
the energy) is wider and the product of the potentials 
is larger than in the other region, therefore contributing 
only a subleading exponential behavior. For this rea- 
son we will only be concerned with the region between 
x*p, = -l/4-fc;/(27r) and x* = -3/4 fc;/(27r) and the 
overlap integral between the two wavefunctions is 



exp 



' exp 



-an 



dx ' 



-2ttx — —TT + k* 
2 " 



-an 1 - 4 , 
16 \ KAao 



2/3^ 



-2tTX TT 

2 



(22) 



The first minimum away from the F point along the 
FY line develops for an > 7. For larger values of a a 
more minima in the band structure can be found (for ex- 
ample, see Fig. ||), increasing the density of states close 
to zero energy, and those band gaps can be computed 
in an analogous fashion to find more exponentially small 
energy splittings. 

Another feature of the band structure that we can un- 



derstand with this simple model is the width of the lowest 
positive and negative energy bands along the FX direc- 
tion. Once again, to find a variational estimate of this 
band width, we need to compute the overlap integral be- 
tween two low-energy wavefunctions. This time, we can 
start from the two zero-energy states at the F point and 
calculate the FX band assuming that it is a narrow tight- 
binding band arising from the overlap of these two states. 
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Again, for large au, we expect the width of this band to 
be exponentiaUy close to zero. The lowest positive and 
negative energy bands along the FX direction in the Bril- 
louin zone can be calculated in the tight-binding approx- 
imation 

£{k,)=±AErxsm^ (23) 

where AE-px is given by the overlap integral of the two 
zero-energy wavef unctions. For large ajj, we can again 
use the WKB approximation to find the asymptotic lead- 
ing behavior of the bandwidth 



showed in the previous section, the potentials U^{x) and 
U^{x) in the A ABB gauge are identical (see eq. (p^), 
therefore the triangular potential wells corresponding to 
a line of A or i? vortices are not spatially separated, as in 
the ABAB case, and the wavefunctions localized around 
the two vortices always have a finite overlap. The nature 
of the asymptotic limit for large a_D is then very differ- 
ent from the case discussed above. We do not expect 
to find an exponential dependence of the features of the 
energy spectrum on the anisotropy, in either one- or two- 
dimensions. Numerical calculations confirm this claim, 
as will be shown in the following section. 



IV. NUMERICAL RESULTS 




FIG. 5. Band structure of the one-dimensional model 
in the ABAB gauge for a square lattice with anisotropy 
ctD ~ vf/va ~ 8, 15, 35 along the FY direction (solid, dashed 
and dot-dashed curve respectively). Only positive energy 
bands are plotted for clarity, negative energy bands can be ob- 
tained through particle-hole symmetry.Energies are in units 
of El — hvp/d. 

To summarize the asymptotic analysis of the one- 
dimensional model in the vortex lattice geometry per- 
formed in this section, we have explained the parallel 
energy bands in the FY direction close to the F point 
and their slope. We have also calculated the leading ex- 
ponential dependence on ao of the FX bandwidth and 
of the minima of the FY bands close to zero energy. In 
the next section, we will describe the numerical solution 
of the onc-dimcnsional model and compare these ana- 
lytic results to the numerical data. These results are im- 
portant to understand the right crossover scale between 
the quantum mechanical and semiclassical regimes of the 
density of states, as discussed in the introduction and in 
our previous publicationQ. 

Very different results would be expected in a singular 
gauge which does not spatially separate the low-energy 
wavefunctions, for example the AABB gauge. As we 



A. Bravais vortex lattice: ABAB gauge 

We have studied the one-dimensional Hamiltonian ( [lO| ) 
numerically discretizing it on a real-space grid in one 
dimension. The reduced dimensionality compared to 
the Bogoliubov-de Gennes operator (H) allows for much 
smaller scale numerical calculations and is the key to ac- 
cessing the large- anisotropy regime. We choose a real- 
space representation of the one-dimensional model in or- 
der to be able to take advantage of sparse matrix diago- 
nalization algorithms (see Appendix B for further details 
on the numerical methods), although we checked some of 
our results against a diagonalization of the Hamiltonian 
( |To| ) in momentum space without finding aay appreciable 
difference. Knapp, Kallin and Berlinskju diagonalized 
the one-dimensional model in the ABAB gauge in mo- 
mentum space and did not find any appreciable difference 
from our results in the range of overlap. 

Let us start discussing the numerical diagonalization 
in the case of the ABAB singular gauge. From the an- 
alytical results discussed in the previous section and our 
previous numerical investigation of the two-dimensional 
Bogoliubov-de Gennes equationQ, we expect to find 
points in the Brillouin zone where the energy spectrum 
has an eigenvalue which is exponentially close to zero. 
The simplicity of the numerical implementation of the 
one-dimensional model allows us to resolve these eigen- 
values in detail and to compare them to our previous 
results in two dimensions. 

Scanning the vortex lattice Brillouin zone we have 
computed the band structure for several values of the 
anisotropy ratio a d ■ The quasiparticle bands in a square 
vortex lattice are symmetric under the exchange of kx — > 
—kx or ky — > —ky, hence only positive k^ or ky have 
generally been considered. Also, particle-hole symmetry 
holds at each point in the Brillouin zone, therefore only 
positive energies need to be taken into account. In the 
one-dimensional model, the potentials U^'^ are periodic 
functions of x and so the band structure is periodic with 
respect to kx^ with period 27r. On the other hand, the 
periodicity with respect to ky is lost, as can be seen, for 
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example, in Fig. ^ where the band structure of the one- 
dimensional model k = (0, 0) (F point) and k = (0, tt) (Y 
point) and the band structure between k = (0, — 27r) and 
fc = (0, —tt) (Y point) are plotted for anisotropy = 12 
respectively with a solid and dashed curve. Note that the 
two branches, differing by a reciprocal lattice vector of 
the two-dimensional vortex lattice Q = (0, —2tt), cross at 
the Y point and at non-symmetry points. These crossings 
become avoided when the non-translational invariance of 
the full two-dimensional Bogoliubov-de Gennes Hamil- 
tonian (|^) is taken into account, but for large enough 
anisotropy, the resulting gaps can be extremely small. 
Also, these crossings and the changes they induce in the 
density of states occur at finite energies and are not too 
crucial in the understanding of the nature of the low- 
energy quasiparticle spectrum. 
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FIG. 6. Band structure of the two-dimensional linearized 
Bogoliubov-de Gennes equation for a square lattice with 
anisotropy ratio ao = vf/va = 8,15,35 (solid, dashed and 
dot-dashed respectively) along the FY direction. Only posi- 
tive energy bands are plotted for clarity and negative energy 
bands can be obtained through particle-hole symmetry. The 
gaps at the F point are fictitious, and are due to the chosen 
lattice discretization scheme (Wilson fermions) of the Hamil- 
tonian. In the an = 15 case, only the two lowest bands are 
plotted, while all the bands with energy E < hvp/d are plot- 
ted for qd = 8 and 12. Energies are in units of Ei — Tivp/d. 

In Fig. 1^ we have plotted the one-dimensional band 
structure for a square vortex lattice with anisotropy ra- 
tio ao — 8,15,35 along the FY axes. This plot should 
be compared to the two-dimensional energy bands calcu- 
lated in the same direction in the Briilouin zone and plot- 
ted in Fig. ^, previously published irfl. The small gaps at 
the F point in Fig. ^ are fictitious and are are due ta,the 
choice of real-space discretization (Wilson fermionsEj) of 
the linearized Bogoliubov-de Gennes equation in two di- 
mensions. In one dimension, we use a different approach 
(staggered fermionstj) which preserves the Dirac node 
at the F point in the discrctized equations, described in 
more detail in Appendix |l|. Similarly to what we ob- 



served for the two-dimensional Hamiltonian (g) , for large 
enough anisotropy, the energy spectrum shows new low- 
energy states around a wavevector close to the Y point. 
As we mentioned before, the degeneracies observed at 
the Y point are split as soon as one introduces the non- 
translational invariance of the Bogoliubov-de Gennes op- 
erator. 

The strong dependence of the minimum (away from 
the F point) of the lowest-energy band along the FY di- 
rection (closest to the Y point) on the anisotropy is 
plotted in Fig. |^, along with the asymptotic expression 
computed in the previous section (the pre-exponential 
factor has been fitted to the one-dimensional numerical 
data) and the same data from a fully two-dimensional 
calculation for values of q;d < 15. The accuracy of the 
asymptotic expression is remarkably good even for rela- 
tively small values of ao and the two-dimensional data 
points are certainly consistent with the one-dimensional 
approximation. 
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FIG. 7. Energy of the minimum of the positive low- 
est-energy band (away from the F point) along FY as a func- 
tion of anisotropy an- The circles (o) represent the results 
of the one-dimensional model while the pluses (-I-) represent 
the two-dimensional one. The solid curve is the asymptotic 
expression calculated in Sect. Ill, with the pre-exponential 
factor fitted to the one-dimensional numerical data. Energies 
are in units of E\ = hvp/d. 

We have also computed the band structure in the FX 
direction and studied its dependence on the anisotropy 
ratio ao- Even for relatively modest values of the 
anisotropy, the lowest energy band is very well approxi- 
mated by the tight-binding expression 

£ik,) = AErxsin^, (25) 

as outlined in the previous section. In Fig. |^, the band 
structure of the one-dimensional model is plotted for 
anisotropy = 4, and it is indistinguishable from the 
tight-binding expression above. We have computed nu- 
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merically AErx over a wide range oi an and the results 
are shown in Fig. ^. The sohd curve is the exponential 
behavior computed in Sect. III. The asymptotic result 



agrees well with the numerical data for large an- 




r X 
FIG. 8. Band structure for the one- dimensional model with 
anisotropy an = vf /v/ Delta = 4 along the FX direction. 
This curve is indistinguishable from the tight-binding expres- 
sion E = 0.15 sin(fci,/2), plotted on the same scale. Energies 
are in units of E-i — hvp/d. 
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FIG. 9. Width of the lowest-energy band along FX as 
a function of anisotropy an- The solid curve represents 
the asymptotic expression calculated in Eq. (p^, with the 
pre-exponential factor fitted to the one-dimensional numeri- 
cal data in the ABAB gauge. The circles (o) are data points 
obtained by diagonalizing numerically the one-dimensional 
model in the ABAB gauge, while the pluses show the 
bandwidth in the FX direction of the two-dimensional Hamil- 
tonian also in the ABAB gauge, for moderate anisotropy. 
The last two data sets plotted show the same energy band- 
width for the one- and two-dimensional model in the AABB 
gauge, plotted with squares (□) and crosses (x), respectively. 
Energies are in units of hvp/d. 



B. Bravais vortex lattice: AABB gauge 



Let us turn to the numerical diagonalization of the 
Bogoliubov-de Gennes Hamiltonian for a square vortex 
lattice in the AABB singular gauge. Although it is still 
a Bravais lattice, the unit cell includes four vortices, as 
shown in Fig. ^ and it is a 1x2 rectangle, in dimensionless 
units. The coordinates of the high-symmetry points in 
the first Brillouin zone arc therefore k = (0, 0) (F point), 
k — (tt, 0) (X point) and fe = (0,7r/2) (L point). As we 
discussed in the previous section, we do not expect to 
find any exponentially small features in the energy band 
structures for large values of the anisotropy ratio in 
this singular gauge. 

The different nature of the low-energy states in the two 
singular gauges can be observed in Fig. ^ where we plot 
the bandwidth of the lowest FX energy band {AErx, as 
defined in equation ( p^ ) in the previous section). No- 
tice the qualitatively different behavior in the two singu- 
lar gauges we considered. While the bandwidth in the 
ABAB gauge follows very closely the asymptotic expo- 
nential dependence on the anisotropy ao, the same en- 
ergy scale in the AABB gauge decays in a much slower 
fashion and does not seem to reach an exponential decay 
in the anisotropy range considered. In both cases, the 
one- and two-dimensional calculations agree reasonably 
well with each other. 

We have studied both the one- and two-dimensional 
Hamiltonians in the AABB gauge. Even though the en- 
ergy bands cannot be computed (even just in the asymp- 
totic limit of large anisotropy) by a simple approach as in 
the ABAB case, the one-dimensional model gives some 
information about the more involved two-dimensional 
Hamiltonian. In particular, the one-dimensional calcu- 
lation of energy bands in the region of the Brillouin zone 
close to the F point is in reasonably agreement with 
the two-dimensional band structure. More generally, the 
translational invariance in the direction of smearing of 
the vortices implies that the different branches of the 
bands in the FL direction — differing by a reciprocal 
lattice vector Q = (0, — tt) — can cross at the L point 
and at non-symmetry points. These crossings, just like 
in the ABAB gauge, become avoided when the trans- 
lational invariance is broken, as in the two-dimensional 
case. Unlike in the ABAB gauge, the magnitude of 
the resulting energy gaps does not become exponentially 
small as we increase the anisotropy ratio. This makes 
the one-dimensional model not as useful for quantita- 
tive calculations of properties of the full two-dimensional 
Bogoliubov-de Gennes Hamiltonian as in the ABAB sin- 
gular gauge. 

The one- and two-dimensional band structures are 
plotted in Figs. |l^ and ^ The energy bands are com- 
puted in the FX and FL directions for anisotropics ao — 
12 and 20. Although the one-dimensional model fails 
to reproduce all the details of the two-dimensional band 
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structure away from the T point even for relatively large 
values of the anisotropy, the one-dimensional bands allow 
to sketch the qualitative behavior of the two-dimensional 
ones, once the band crossings become avoided. Finally, 
notice that, while in the ABAB gauge there was a min- 
imum developing in the lowest energy band already for 
a_D — 8, in the A ABB gauge the anisotropy ratio has 
to be much larger {ao = 20) before we start finding an 
inflection in the lowest energy band for both the one- and 
two-dimensional Bogoliubov-de Gennes operators. 
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FIG. 10. Band structure for a square vortex lattice in 
the AABB singular gauge with an = vf/va ~ 12. The 
solid line is the spectrum of the two-dimensional Bogoli- 
ubov-de Gennes equation, while the dashed line represents 
the one-dimensional model. Only positive energy bands are 
plotted for clarity, negative energy ones can be obtained by 
particle hole symmetry. Energies are in units of Ei = Tivp/d. 

In conclusion, even though the band structures in the 
ABAB and AABB gauges are qualitatively similar close 
to the center of the Brillouin zone, they differ signifi- 
cantly near the edges of the Brillouin zone. While in the 
ABAB case all relevant energy scales depend exponen- 
tially on the anisotropy ratio ud, the same is not true 
for the AABB gauge. This behavior is well understood 
in terms of the simplified one-dimensional model in both 
gauges, for large values of the anisotropy ratio au- 



but let us consider, for concreteness, the case where the 
distance between nearest A and B vortices is (2\/2/5)d, 
which corresponds to vortex coordinates in the unit cell 
of (±1/5, ±1/5) in dimensionless units. We can use the 
one-dimensional potentials computed in equations jl^ ) 
and jl^ ) with = 1/5 to determine the one-dimensional 
band structure. The energy bands along the -Y to Y di- 
rection for anisotropy ratios ao = 1 and 8 are plotted 
in Figs. 12 and |l^, respectively. 

-kr and h 



First of all, we notice 
fc„ symmetries of the 



that the k 

energy spectrum are broken. Fig. O exemplifies how the 
one-dimensional very effectively captures the qualitative 
features of the two-dimensional band structure also for a 
non-Bravais lattice in the large anisotropy limit for the 
same reasons discussed above in the Bravais lattice case. 




X r L 

FIG. 11. Band structure for a square vortex lattice in 
the AABB singular gauge with an = vf/va = 20. The 
solid line is the spectrum of the two-dimensional Bogoli- 
ubov-de Gennes equation, while the dashed line represents 
the one-dimensional model. Only positive energy bands are 
plotted for clarity, negative energy ones can be obtained by 
particle hole symmetry. Energies are in units of Ei = Tivf/cL. 



C. Non-Bravais vortex lattice 

The Bravais lattice nature of the vortex lattices we 
considered so far ensures the gaplessness of the spectrum 
at the center of the Bruillouin zoneQ. In the ABAB sin- 
gular gauge, we can consider a unit cell with two vortices 
again, but relax the constraint that they sit equidistantly 
from each other, along the diagonal — i.e. the magni- 
tude of iio, as defined in Fig. |l[ does not need to be 
equal to •\/2rf/4. Let us assume that the two sublat- 
tices A and B are still square lattices with spacing d, 



We want to emphasize here that, unlike for a Bravais 
lattice of vortices, the energy spectrum here is gapped. 
These gaps where discussed in a perturbation theory 
framework inQ. What is interesting to note is that the 
evaluation of these gaps in the one-dimensional model 
is very consistent with the two-dimensional calculation, 
even though the handling of the discretization of the lin- 
earized Bogoliubov-de Gennes equation is different in the 
one- and two-dimensional case, as is explained in more 
detail in Appendix B. 
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FIG. 12. Band structure for a non-Bravais square vor- 



tex lattice in the ABAB singular gauge from -Y to Y with 
an = vf/va ~ 1 and vortex coordinates in the unit 
cell (±1/5, ±1/5). The solid line is the spectrum of the 
two-dimensional Bogoliubov-de Gennes equation, while the 
dashed line represents the one-dimensional model. Only posi- 
tive energy bands are plotted for clarity, negative energy ones 
can be obtained by particle hole symmetry. Energies are in 
units of El — hvp/d. 
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FIG. 13. Band structure for a non-Bravais square vor- 
tex lattice in the ABAB singular gauge from -Y to Y with 
Old = Vf/va ~ 8 and vortex coordinates in the unit 
cell (±1/5, ±1/5). The solid line is the spectrum of the 
two-dimensional Bogoliubov-de Gennes equation, while the 
dashed line represents the one-dimensional model. Only posi- 
tive energy bands are plotted for clarity, negative energy ones 
can be obtained by particle hole symmetry. Energies are in 
units of El — hvF/d. 



V. CONCLUSIONS 

We have studied the large Fermi-velocity anisotropy 
ao = vf/va hmit of the quasiparticle spectrum of a 
d— wave superconductor in the mixed state. The vortex 
geometry we considered is that of a square vortex lattice 
rotated by 45° from the anisotropy axes. The original 
linearized Bogoliubov-de Gennes equation was mapped 
onto a Dirac Hamihonian in an effective periodic vec- 
tor and scalar potential using the Franz-Tesanovic gauge 
transformation. We have studied the dependence of the 
energy spectrum on the choice of such singular gauge 
transformation. Physically, one might expect the gauge 
choice not to matter. However, although the Bogoliubov- 
de Gennes operator considered in this paper is indeed 
gauge invariant away from the vortices, the singularity 
at the vortex site introduces boundary conditions that 
break the invariance of the energy spectrum on the choice 
of singular gauge. 

In the large anisotropy limit, we have solved a one- 
dimensional variant of the Franz-Tesanovic model ob- 
tained by smearing the vortices in the "hard" direction. 
We computed the one-dimensional potentials associated 
with two different choices of singular gauges in the case of 
a Bravais square vortex lattice. In one case (the ABAB 
gauge), we found that the low-energy localized states 
around each vortex type {A or B) in the one-dimensional 
potential wells are physically separated. The model can 
therefore be solved in the WKB approximation and sev- 
eral features of the energy band structure are shown to 
depend exponentially on the anisotropy ratio ud- We 
compared the computed energy spectrum to the results 
of a numerical diagonalization of the Hamiltonian ([ic|). 
In both cases we found that both the minima of the low- 
est energy-band in the FY direction and the bandwidth 
of the lowest energy-band in the FX direction decay ex- 
ponentially to zero as a function of the anisotropy au- 
The asymptotic calculation agrees very well with the nu- 
merical results. 

Different results are obtained in the AABB gauge. In 
this case, the low-energy states localized around the vor- 
tices of different types are not spatially separated. This 
results in a finite overlap integral for any value of the 
anisotropy and, therefore, no exponentially small features 
in the band structure. The one- and two-dimensional 
numerical spectra in this gauge are not as similar as in 
the ABAB gauge, but one can still extract information 
like degeneracies or other global properties of the two- 
dimensional band structure from the one-dimensional 
calculations. Of course, the detailed spectra will not be 
equal even in the asymptotic large ao limit. 

Finally, we have distorted the square vortex lattice 
in the ABAB gauge into a non-Bravais vortex lattice. 
Again, we find good agreement between the one- and 
two-dimensional calculations of the energy spectra. In 
particular, we find a gap in the lowest energy band close 
to the F point (exactly at the center of the Brillouin zone 
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only in the ao = 1 case, otherwise the —k^ and 

ky —ky symmetry of the energy bands is broken) in 
both numerical calculations. 
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APPENDIX A: WKB APPROXIMATION FOR 
THE ONE-DIMENSIONAL MODEL 

We want to find the eigenfunctions for a given energy 
E of the one-dimensional model in the WKB approxima- 
tion and derive the appropriate Bohr-Sommerfeld quan- 
tization condition. Starting from the one-dimensional 
Hamiltonian (]To|), let us define the potentials ^^{x) = 
ky + U^{x) and ^^{x) = -ky + U^{x) to simplify nota- 
tions. If we fix the energy we need to find a solution 
to the system of differential equations 



^^{x)u{x) + (^a, + ^fc,) v{x) = Eu{x) 
[ih^- + -k^^) + '^''{^)v{x) = Ev{x). 



(Al) 



The derivation of the WKB form of the wavrfunctions 
proceeds along standard lines (see for examplec3). Sub- 
stituting the expansion 



U[x) — e ^n=o ' o 

V[X) = e ^n=0 ' o 



(A2) 



in the previous system of equations, to find the asymp- 
totic behavior of the wavefunctions for large ao- Ex- 
pand ing and equating order by order in q;d in equation 
(Al), we can calculate the first two terms in the series 
(A2), namely 5'o,i(a;) and Tq^i{x). These are the leading 
terms in the expansion for large a d and the truncation of 
the series (A2) to this order is the WKB approximation. 
In particular, Sq{x) and Tq{x) turn out to be 

^0(2:) = To{x) - ± / dx' ^(i;-$^(a;'))(-B-$^(x')) 



(A3) 

which gives the leading exponential behavior and is the 
key expression for the Bohr-Sommerfeld quantization for- 
mula. Every other term in the series expansion leads 
to algebraic corrections, which we have not taken into 
account in this work. For completeness, we state the 
expression of the WKB wavefunctions, including the al- 
gebraic prefactor: 



u{x) 



1/4 



,±^k^x ±^»D P dx' ^(B-*-* (2;')) (2;')) 



(4) 



Obviously, these oscillating expressions are correct in 
the classically allowed regions, while the exponential fac- 
tors become real in the classically forbidden regions. 

APPENDIX II: NUMERICAL 
IMPLEMENTATION OF THE 
ONE-DIMENSIONAL MODEL 

To diagonalize the Hamiltonian ( |lO| ) we have used 
a real-space representation to take advantage of fast 
sparse-matrix diagonalization algorithms and to con- 
trol the smoothing of vortex cores over a finite region 
of space, which cannot be avoided in a momentum- 
space representation. The main drawback of discretiz- 
ing a Dirac-type Hamiltonian in real space is the, ap- 
pearance of the "Fermion doubling problem" (seeEj for 
a thorough discussion of possible ways of discretizing a 
Dirac Hamiltonian). The real-space representation of a 
Dirac Hamiltonian leads to the introduction of 2^ — 1 
spurious large- wavevector, low-energy modes (additional 
Fermions, hence the name), where D is the dimension- 
ality of the system. In our case, the system being one 
dimensional, there is only one spurious mode, which we 



can eliminate using the staggered-fermion approachEH. If 
we discretize the vortex lattice unit cell with a mesh of 
size h = 1/ (N — 1) and index the resulting mesh with an 
integer n ranging from —N/2 to N/2 so that x — n/N, 
we can place a single-component Fermi field ^(n) on each 
site. To identify a single two-component Dirac field, we 
can decompose the mesh into an even and an odd sub- 
lattice and define 

l/(n) = ^(n), n even , . 

W{n) = ^{n), n odd. ^ ' 

The staggered fermion method avoids the Fermion dou- 
bling problem by doubling the size of the unit cell in real 
space, thus halving the size of the Brillouin zone (and 
overall number of degrees of freedom). Note that the 
relevant unit cell here is unrelated to the vortex lattice 
unit cell, it is simply the unit element of the mesh with 
which we discretize the continuum equations. Defining 
the function 

r(r,\ - f U^{n) + ky, n even , , 

^^'''-\U^{n)-ky, nodd 

the discretized eigenvalue problem takes the form 
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1 + 



ao 2 



(3) 



and the correct continuum limit is achieved by let- 
ting h ~* 0, without any further caveats. An importapt 
advantage of staggered ferrwons over Wilson fermionsE2l 
(which we used previousljfl) is that Wilson approach 
breaks the chiral symmetry of the translationally invari- 
ant Hamiltonian, introducing a fc— dependent mass term 
to lift the spurious modes at the edges of the mesh Bril- 
louin zone to high energy. By explicitly breaking the 
chiral symmetry, fictitious gaps can and do appear in 
the energy spectrum at the center of the vortex lattice 
Brillouin zone thus making it difficult to perform a pre- 
cise study of almost dispersionless bands, which is one of 
the goals of the present work. In principle, these gaps 
could be eliminated to arbitrary precision by carefully 
tuning suitable counterterms but the procedure (having 
to be performed numerically) is slow and not very accu- 
rate. The staggered fermion approach bypasses all these 
difficulties by preserving explicitly the chiral symmetry 
of the one-dimensional Hamiltonian in the lattice rep- 
resentation, if there was one initially in the continuum 
representation. 
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